Fragestellung

Wie kann die Mehrebenenstruktur der Sleep Study Datensatzes beschrieben werden?

Lösungshinweis 1

Mit head(sleepstudy) erkennt man, dass der Datensatz nur aus AV, UV, Subjectidentifier besteht. Also muss können maximal zwei Ebenen im Datensatz vorkommen.
Bildet man ein Tabelle mit Days und Subject kann man entscheiden ob eine perfect iode rimperfecte Hierarchie vorliegt.

Lösungshinweis 2

Eine solche Tabelle erhält man mit table(sleepstudy$Days, sleepstudy$Subject)

Lösung

Es handelt um eine perfekte Zweiebenenstruktur mit folgenden Nestgrößen und folgender Nestanzahl:

library(lme4) #enthält den Datensatz
# Nestanzahl
length(unique(sleepstudy$Subject))

## [1] 18

# Nestgrößen
sleepstudy %>% 
  group_by(Subject) %>% 
  summarize(Nestgröße = n())

Fragestellung

Der ICC hilft dabei empirisch zu erkennen wie stark sich die Mehrebenenstruktur in der abhängigen Variable ausdrückt. Dazu schätzt man ein Intercept-Only Modell und setzt die Intercept-Varianz ins Verhältnis zur Gesamtvarianz.

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]}, \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \]

Lösungshinweis 1

Die Syntax für das Intercept-Only Modell lautet lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)

Lösungshinweis 2

Richtig schnell geht es mit der Funktion icc() aus dem Package {performance}.

Lösung

library(performance)
library(lme4)
mod01_sleepstudy <- lmer(Reaction ~ 1 + (1|Subject), data = sleepstudy)
icc(mod01_sleepstudy)

## # Intraclass Correlation Coefficient
## 
##      Adjusted ICC: 0.395
##   Conditional ICC: 0.395

Dieser ICC ist vglw. groß! Multi-Level Modelling ist also nicht nur theoretisch impliziert sonder auch empirisch.

Fragestellung

Um später entscheiden zu können inwiefern Random Interept und Slope empirisch gerechtfertigt sind, kann man zunächst ein Single Level Model schätzen um dieses später mit den komplexeren Multi-Level Modellen zu schätzen.

Lösung

mod02_sleepstudy <- lm(Reaction ~ Days, data = sleepstudy)
summary(mod02_sleepstudy)

## 
## Call:
## lm(formula = Reaction ~ Days, data = sleepstudy)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -110.848  -27.483    1.546   26.142  139.953 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  251.405      6.610  38.033  < 2e-16 ***
## Days          10.467      1.238   8.454 9.89e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 47.71 on 178 degrees of freedom
## Multiple R-squared:  0.2865, Adjusted R-squared:  0.2825 
## F-statistic: 71.46 on 1 and 178 DF,  p-value: 9.894e-15

Das entspricht einem Modell \[ \operatorname{Reaction} = {\color{#8cd000}{\beta_{0}}} + {\color{#8cd000}{\beta}}_{1}(\operatorname{Days}) + {\color{#8cd000}{\epsilon}} \] mit den Koeffizienten \[ \operatorname{\widehat{Reaction}} = 251.41 + 10.47(\operatorname{Days}) \]

Fragestellung

Wie ändern sich die Koeffizienten des Modells, wenn man zulässt, dass Intercepts zwischen den Clustern variieren?

Lösungshinweis

Das Intercept kann ja als konstanter Prädiktor mit dem Wert 1 aufgefasst werden. Daher lautet die Syntax für das Random Intercept (1|Subject) (lies: »das Intercept darf über die Clustervariable Subject variieren«).

Lösung

mod03_sleepstudy <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy)
summary(mod03_sleepstudy)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

Das entspricht einem Modell

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(\mu_{\alpha_{j}}, \sigma^2_{\alpha_{j}} \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] mit den Koeffizienten

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1}}(\operatorname{Days}), \sigma^2 \right) \\ \alpha_{j} &\sim N \left(0, 37.12 \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] Die Modelle vergleichen kann man mit der Syntax

tab_model(mod01_sleepstudy, mod02_sleepstudy, mod03_sleepstudy)

Die Modelle gegeneinander Testen \(H_0: \text{beide Modelle klären gleich viel Varianz auf}\) kann man mit dem anova()-Befehl:

anova(mod03_sleepstudy, mod02_sleepstudy)

## refitting model(s) with ML (instead of REML)

Unterschied ist signifikant, also steigert das Random Intercept die Modellpassung überzufällig.

Fragestellung

Wie ändern sich die Koeffizienten des Modells, wenn man zulässt, dass Intercept und Slopes variieren?

Lösungshinweis

Die Syntax für ddie Random Effects lautet (1 + Dayes|Subject) (lies: »das Intercept und die Slopes dürfen über die Clustervariable Subject variieren«).

Lösung

mod04_sleepstudy <- lmer(Reaction ~ Days + (1 + Days|Subject), data = sleepstudy)
summary(mod03_sleepstudy)

## Linear mixed model fit by REML ['lmerMod']
## Formula: Reaction ~ Days + (1 | Subject)
##    Data: sleepstudy
## 
## REML criterion at convergence: 1786.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2257 -0.5529  0.0109  0.5188  4.2506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Subject  (Intercept) 1378.2   37.12   
##  Residual              960.5   30.99   
## Number of obs: 180, groups:  Subject, 18
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 251.4051     9.7467   25.79
## Days         10.4673     0.8042   13.02
## 
## Correlation of Fixed Effects:
##      (Intr)
## Days -0.371

Das entspricht einem Modell

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{Reaction}_{i} &\sim N \left(\alpha_{j[i]} + \beta_{1j[i]}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &\mu_{\alpha_{j}} \\ &\mu_{\beta_{1j}} \end{aligned} \end{array} \right) , \left( \begin{array}{cc} \sigma^2_{\alpha_{j}} & \rho_{\alpha_{j}\beta_{1j}} \\ \rho_{\beta_{1j}\alpha_{j}} & \sigma^2_{\beta_{1j}} \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] mit den Koeffizienten

## Warning: Colorization of greek notation not currently implemented for merMod
## models

\[ \begin{aligned} \operatorname{\widehat{Reaction}}_{i} &\sim N \left(251.41_{\alpha_{j[i]}} + 10.47_{\beta_{1j[i]}}(\operatorname{Days}), \sigma^2 \right) \\ \left( \begin{array}{c} \begin{aligned} &\alpha_{j} \\ &\beta_{1j} \end{aligned} \end{array} \right) &\sim N \left( \left( \begin{array}{c} \begin{aligned} &0 \\ &0 \end{aligned} \end{array} \right) , \left( \begin{array}{cc} 24.74 & 0.07 \\ 0.07 & 5.92 \end{array} \right) \right) \text{, for Subject j = 1,} \dots \text{,J} \end{aligned} \] Die Modelle vergleicht man mit der Syntax

tab_model(mod01_sleepstudy, mod02_sleepstudy, mod03_sleepstudy, mod04_sleepstudy)

Die Modelle gegeneinander Testen \(H_0: \text{beide Modelle klären gleich viel Varianz auf}\) kann man mit dem anova()-Befehl:

anova(mod04_sleepstudy, mod03_sleepstudy)

## refitting model(s) with ML (instead of REML)

Unterschied ist signifikant, also steigert der zusätzliche Random Slope die Modellpassung überzufällig.